Synchronization of optically self-assembled nanorotors

Self-assembly of nanoparticles by means of interparticle optical forces provides a compelling approach toward contact-free organization and manipulation of nanoscale entities. However, exploration of the rotational degrees of freedom in this process has remained limited, primarily because of the predominant focus on spherical nanoparticles, for which individual particle orientation cannot be determined. Here, we show that gold nanorods, which self-assemble in water under the influence of circularly polarized light, exhibit synchronized rotational motion at kilohertz frequencies. The synchronization is caused by strong optical interactions and occurs despite the presence of thermal diffusion. Our findings elucidate the intricate dynamics arising from the transfer of photon spin angular momentum to optically bound matter and hold promise for advancing the emerging field of light-driven nanomachinery.


Supplementary Text Optical forces calculation
The calculation of the optical forces acting on the nanoparticles is a basic ingredient in our simulation of the nanoparticle motion.The electromagnetic fields have been calculated using the Ansys Lumerical FDTD program as well as using T-matrix methods combined with extended Mie theory (48,49).
Once we know the electric and magnetic fields  ⃗ and  ⃗ ⃗ on a surface  enclosing a particle, we can calculate the force  on the particle and the torque  ⃗⃗ about its center,  0 , in terms of surface integrals involving the Maxwell stress tensor  ⃡ through (here and in the following ℜ denotes the real part) and The integration, and hence the point  and the corresponding outward-pointing normal vector  ̂, runs over a surface enclosing the particle, see for example (50) for further details.The elements of the stress tensor are given by where the indices  and  denote the Cartesian coordinate directions and   is the relative permittivity of the medium surrounding the particle (water in our case).We assume that the medium surrounding the particles is non-magnetic.

Simulation of the particle dynamics
With the knowledge of the optical forces we use them in a simulation of the particle motion.In these calculations we model the particles as prolate spheroids with a semiaxis of length  along its cylindrical symmetry axis and a semiaxis  in the perpendicular direction.Their translational motion is restricted to a plane while they can also rotate around one axis that is orthogonal to that plane and to the symmetry axis of the particle.The translational motion follows Newton's second law, where, however, the first inertial term, involving the particle mass   can be neglected.
The second term on the left hand side describes a viscous friction force that varies linearly with the velocity   of particle  (given that the flow around the particle is associated with a very low Reynolds number).The frictional forces are calculated using the formulas obtained by Perrin (51) a long time ago to describe the viscous friction of spheroids and ellipsoids in a liquid.On the right hand side  opt () is the optical force, depending on the positions { } = ( 1 ,  2 , … ) and orientations { ⃗ } = ( ⃗ 1 ,  ⃗ 2 , … ) of all the particles.The last term on the right hand side of Eq. ( 4) yields a stochastic force due to thermal fluctuations.
In addition to causing friction and giving rise to stochastic forces, the water surrounding the particles can mediate hydrodynamic interactions between them.As one particle spins around its own axis, it sets water in motion generating a flow at the other particles.However, the speed of the flow decays as 1  2  ⁄ with increasing distance  from the spinning particle (and the decay is even faster than that in case there is a fixed substrate nearby).
This means that in our case the velocity of the flowing water decreases by at least two orders of magnitude going from the surface of a spinning particle to a neighboring particle.Given the linear relation between torques and spinning frequencies, this also gives us an estimate of the ratio between the torque caused by hydrodynamic interactions and the optical torque.In view of this we conclude that hydrodynamic interactions have very little impact on the rotational motion of the nanoparticles since they only account for torques that are of the order of 1% or less of the optical torques.
In our dynamic simulation then for each time step (typically 2 µs) we calculate a deterministic velocity as a result of the balance between friction and optical forces,   (  ) =  ⃡  −1  opt () , (5) which results in a deterministic displacement during the time-step.To this we add a stochastic displacement (∆ ‖ and ∆ ⊥ ) drawn from a Gaussian distribution set by the temperature and the diffusion coefficients (the diffusion rates along and orthogonal to the particle symmetry axis differ).Then, the new coordinates for a particle at time  +1 =   +  are given by The rotational motion is governed by completely analogous equations with the deterministic angular velocity given by ̇(  ) =  rot −1  opt () .(8) The theoretical heat maps showing the correlations between the scattered light intensity from different particles were calculated using the dynamic simulations based on Mie theory.In these simulations we mimicked the experimental illumination and subsequent scattering of light (figs.S2 and S7) by following the orientation of the particles during time intervals of length   = 100 µs.In the simulation, these time intervals follow each other without intermissions.Given that the scattering of the probing light is primarily caused by dipolar radiation emanating from the long axis dipole of the particles, to a first approximation the scattered electric field scales as   , and the scattered intensity, apart from a constant, is proportional to  2   , see fig.S7, and we will from now on refer to it as normalized intensity.Hence, to calculate the variations in scattering intensity we run the dynamic simulation using a time step of 2 µs.At each time step the orientation of each particle is recorded and then the "flash" average of the normalized intensity is calculated for each particle and each (theoretical) flash.The so obtained heat maps of normalized intensity show a clear correlation between different particles, provided the incident laser power that causes optical binding and drives the rotational motion is sufficiently strong.In simple terms this means that the different nanoparticles tend to have their major axis' oriented in nearby directions, an effect caused by the electromagnetic coupling between them.
To get a quantitative measure of the synchronization of the spinning motion we have calculated the order parameter (here   is the total number of particles) and averaged it in time over the course of the simulations with up to   = 2.5 × 10 6 points in time.The results show values of 〈〉 ranging from 0.81 to 0.91 for the laser powers used in the case of two spinning particles.

Qualitative discussion of the synchronization mechanism
The basic mechanism behind the synchronization of the particles' spinning can be understood within a simple model involving electric dipole-dipole interactions.Similar models have been applied in the analysis of other experiments involving optical interactions and synchronization (24,31).
Consider two point dipoles situated in a plane at a distance  from each other, as illustrated in fig.S6A, with the dimer axis aligned with the x axis.We assume that the dipoles are immersed in water and that typically the distance between them equals the optical wavelength in water  0 / for the incident light driving the dipoles.Assume also, for the moment, that the dipoles are highly anisotropic with a considerable polarizability only in one direction in the plane of the dimer (that is the xy plane).
These dipoles are now forced to spin by an external electric field,  ⃗ 0 =  ̂0 +  ̂0 , (12) equal at both dipoles, and with an implicit time dependence  − .The response of the dipoles are determined by an effective polarizability tensor  ⃡  yielding   =  ⃡   ⃗  , (13) where  denotes the particle number.Note that  ⃗  in the above equation denotes the total field at dipole , also including fields that have been scattered off the other dipole.
Given the highly anisotropic response of the dipole the polarizability tensor depends on the angle   formed between the dipole axis and the dimer axis, hence where  0 is a constant with real and imaginary parts  0 ′ and  0 ′′ respectively, The imaginary part is due to losses, partly through absorption but also, since this is an effective polarizability, as a result of scattering.
The time-averaged (over one period of the electromagnetic oscillation) torque on particle (dipole) 1 can be determined from the simple expression (ℜ and ℑ denote the real and imaginary parts, respectively) which can be written as With one single dipole, the electric fields entering this expression are identical with the incident field  ⃗ 0 .Linearly polarized light will then align the dipole along the direction of the incident field.
For circularly polarized light, on the other hand, the first part of the torque in Eq. ( 17), associated with the real part of the polarizability, vanishes, whereas the second part, proportional to  0 ′′ , yields a torque that turns the particle in the direction of the circular polarization.
Next we look specifically at the dipole-dipole interaction.With a separation between the dipoles equal to the wavelength of light in the medium, this interaction is approximately limited to dipole moments perpendicular to the dimer axis (dipole in the y direction) giving rise to a field in the same direction on the other dipole.That is, the total field at dipole 1 can approximately be written In the particular situation we are considering  is, to a good approximation, real and positive and in the following we treat it as a real and positive constant.We further take the incident field to be  ⃗ 0 =  0 ( ̂+  ̂).Then the correction  ⃗ 1 to the field causes a first-order correction (in ) to the optical torque equal to In the cases we consider here, the real part of the polarizability dominates the imaginary part which means that the leading corrections to the torque caused by particle-particle interactions comes from the first two terms, proportional to  0 ′2 , in Eq. ( 19).The first term in  1 , varying as 2 1 (but still related to interactions as it is proportional to ), tends to align the dipole with the y direction.This is due to the fact that there is an additional electric field in this direction (which is in phase with the y component of the external field).This contribution to the torque will in general cause the dipole to spend more time pointing in the y direction than the x direction while performing its spinning motion.As a result of this phenomenon we see oscillations in the spinning frequencies in the simulations.
The second term in Eq. ( 19), varying as (2 2 − 2 1 ), is even more interesting.We see that if the two angles  1 and  2 differ, the case with a somewhat larger  2 will give an additional positive contribution to the torque  1 that means that  1 will tend to catch up with  2 , and vice versa.This interaction, which is identical in form to the interaction causing synchronization in the Kuramoto model (38) is at the core of the synchronization of the spinning motion of two or more nanoparticles that we see in our experiments and simulations.Moreover, the order parameter defined in Eq. ( 10) is also used as a measure of synchronization in the Kuramoto model.
To make contact with the expressions for the spinning frequencies given in the main text, and the quantities  0 and ∆ used there, we have To check the validity of the above reasoning we have, using the T-matrix formalism, calculated the torques on two gold particles in the form of prolate spheroids with semi-axis of 55 nm and 88 nm, respectively.Figure S6B shows calculated results for the torque in three different scenarios (i) the particle response is limited to the electric dipole along the major axis, (ii) the response is limited to electric dipole response along all three axes, and (iii) we allow for the full response of the particle including all relevant multipoles.For comparison a fourth curve shows the (orientation-independent) torque one would get if the particles were spheres with a radius of 88 nm.
The results show that the assumption of a particle that can only be polarized along one direction is too simple in the sense that it overestimates the torque as well as the corrections to the torque due to particle-particle interactions.Allowing for electric dipole response along all axes, on the other hand, give results in close agreement with those obtained including the full response of the particles.
Figure S6C shows the torque on particle 1 as a function of the orientations of both particles,  1 and  2 .These results, to a large degree, confirm the qualitative predictions of the dipole model above.The torque  1 shows a strong variation with  1 similar to the first term in Eq. ( 19) and also a variation with  2 −  1 similar to the second term in Eq. (19).nanoparticle with semi-axis 55 nm and 88 nm, respectively.The particle is separated from another identical particle by a distance of 810 nm and both particles are illuminated by a plane electromagnetic wave at normal incidence with a wavelength of 1064 nm and intensity 0.1W/(μm 2 ).
Here the other particle has a fixed orientation with its long axis perpendicular to the dimer axis.(C) Color plot of the torque M1 as a function of the orientation of both particles in the dimer.

Movie S2.
Random transitions of four optically bound gold nanorods.

Movie S3.
Random transitions of five optically bound gold nanorods.

Movie S4.
Orbital rotation of two optically bound gold nanorods.

Movie S5.
Orbital rotation of three optically bound gold nanorods.

Movie S6.
Orbital rotation of four optically bound gold nanorods.

Movie S7.
Orbital rotation of five optically bound gold nanorods.Movie S1 to S3 have been slowed down by a factor 100. Movie S4 to S7 have been slowed down by a factor 10.

Fig. S1 .
Fig. S1.Optical properties.Measured ensemble-averaged extinction spectrum of the gold nanorod (Au NR) sample.The red dash line indicates the trapping wavelength of 1064 nm.

Fig. S2 .
Fig. S2.Optical setup.Schematic of the optical setup used to observe self-organization and synchronization of NR optical matter arrays.The NRs are trapped in two dimensions against the upper cover glass of thin a liquid sample cell by a 1064 nm wavelength laser beam.The beam size was truncated by an iris before being focused on the back-focal plane of the objective in order to minimize beam divergence at the sample plane.The polarization of the beam was adjusted to left circular by manually tuning a half-wave and a quarter-wave plate.The NRs are illuminated through an oil-immersion dark-field (DF) condenser with white light from the microscope halogen lamp or red light from a LED.For stroboscopic observation of NR rotation behavior, a linear polarizer is inserted in front of the CMOS camera and the on/off time of the LED is set to 0.1/1.58ms while the camera frame rate is set to 633 Hz.

Fig. S3 .
Fig. S3.Trapping laser focus.(A) Optical image of the laser spot captured by placing a mirror at the trapping plane.(B) Measured (red line) and Gaussian fit (black dashed line) of the intensity profile in the trapping plane.The FWHM of the profile is ~2.5 μm.

Fig. S4 .
Fig. S4.Orbital rotation characteristics.(A) Evolution of NR angle relative to x-axis for arrays of 2-6 close-packed NRs.(B) Orbital frequency versus number of NRs from the data in (A).The incident laser power was ~430 mW.

Fig. S5 .
Fig. S5.Model NRs.(A) NR dimensions used in simulations.FDTD simulations were based on modelling a NR as a cylinder with ellipsoidal caps with Rxy = 50 nm, Rz = 26 nm, and Lz = 166 nm (blue).Mie calculations and Brownian diffusion simulations were based on modelling a NR as a prolate spheroid with minor axis rxy = 55 nm and major axis rz = 88 nm (green).(B) Simulated extinction spectrum of a NR using FDTD.(C) Calculated orientation-averaged extinction spectrum of NRs using Mie theory.

Fig. S6 .
Fig. S6.Variation in optical torque for a dimer of prolate spheroids.(A) Schematic illustration of the particles forming the dimer and their orientations.(B) Calculated torque M1 on a prolate spheroidal nanoparticle with semi-axis 55 nm and 88 nm, respectively.The particle is separated from another identical particle by a distance of 810 nm and both particles are illuminated by a plane electromagnetic wave at normal incidence with a wavelength of 1064 nm and intensity 0.1W/(μm 2 ).

Fig. S7 .
Fig. S7.Methodology to determine NRs' orientations.(A) Schematics of the experiment.The sample is illuminated through a DF oil condenser with NA=1.2-1.43 using unpolarized red light (λ=740 nm) from an LED.The refractive indexes of the used oil, cover glass and water are 1.515, 1.52 and 1.33, respectively.Light scattered in the forward direction pass through a linear polarizer oriented along the x-axis and is collected by a CMOS camera.(B) Simulated DF intensity versus NR orientation.The NR is illuminated with an x-polarized plane wave of wavelength λ=740 nm.The power scattered into the forward angular range between 64.5 and 90° light is analyzed as a function of angle  between the NR long axis and the x-axis.Due to electromagnetic reciprocity, this set up is a good representation of the experiment.(C) Simulated DF intensity versus NR angle .

Fig. S8 .
Fig. S8.Laser power-dependent scattering intensity distributions.(A) Measured distribution of single particle scattering intensities recorded using stroboscopic illumination of optically bound NR dimers trapped using ~100, 320, 430, and 530 mW laser power, respectively.(B) Calculated distribution of single particle scattering intensities by running a dynamic simulation with two optically bound prolate spheroids using 400, 700, 800, and 900 mW laser power, respectively.

Fig. S9 .
Fig. S9.Intensity statistics for dimers of optically bound Au nanospheres.Scatter plots of intensities I1 vs. I2, extracted from >10000 separate images of two pairs of nominally spherical Au nanoparticles and binned into heat-maps, for two different laser powers.The distributions are narrow and almost circular.

Fig. S10 .
Fig. S10.Simulated intensity distribution vs. flash time.The figure shows the cumulated intensity probability distribution perpendicular to the diagonal I1=I2, i.e. across the heat map ridge, from stochastic simulations of an optically bound dimer versus simulated integration times ("flashtime").The shortest flash-time produces a distribution with a sharp cusp at  = 0 .This is consistent with a Boltzmann distribution in angular space ( 1 ,  2 ) based on the nanorod interaction.

Fig. S11 .
Fig. S11.Simulated intensity distribution and randomization tests.The figure shows simulated intensity correlations I1 vs. I2 from stochastic simulations of an optically bound NR dimer.(A) Original data; (B) Result when I1 is reduced by a factor 0.6; (C) Result when I1 is reduced by a factor 0.6 and the particle identities are scrambled; (D) Result when I1 is reduced by a factor 0.6, the I2 data is delayed by 1 s to cancel all correlations, and the particle identities are scrambled.The data shows that the positive co-variance is only marginally affected by a moderate variation in particle scattering amplitude and by scrambling the particle identities (as in the experiments).

Fig. S12 .
Fig. S12.Statistical analyses of two additional NR trimers.(A) DF microscopy using stroboscopic redlight illumination combined with polarization-filtered imaging reveal intensity fluctuations due to variations in NR orientation with the respect to the polarizer.Snapshots of a trimer using stroboscopic red-light illumination combined with polarization-filtered imaging together with schematics indicating possible NR orientations compatible with the recorded particle intensities.(B, C) Grid heat maps of the intensity correlations for two additional trimers measured using an applied laser power of ~430 mW.

Fig. S13 .
Fig. S13.Statistical analyses of two additional NR tetramers.(A, B) Grid heat maps of the intensity correlations for two additional trimers measured using an applied laser power of ~430 mW.